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ABSTRACT 

Parameter studies are conducted using the Euler and potential flow equation models for steady and unsteady flows in 
both two and three dimensions. The Euler code is an implicit, upwind, finite volume code which uses the Van Leer method 
of flux- vector* split ting which has been recently extended for use on dynamic meshes and maintains all the properties of 
the original splitting. The potential flow code is an implicit, finite difference method for solving the transonic small- 
disturbance (TSD) equations and incorporates both entropy and vorticity corrections into the solution procedure thereby 
extending its applicability into regimes where shock strength normally precludes its use. Parameter studies resulting in 
benchmark type calculations include the effects of spatial and temporal refinement, spatial order of accuracy, far field 
boundary conditions for steady flow, frequency of oscillation, and the use of subiterations at each time step to reduce 
linearization and factorization errors. Comparisons between Euler and potential flow results are made as well as with 
experimental data where available. 

INTRODUCTION 

Considerable progress has been made over the past two decades on developing computational fluid dynamics (CFD) 
methods for steady and unsteady flow analysis. Two recent surveys which summarize much of this progress are given in 
Refs. 1 and 2. Jameson, 1 for example, has discussed successes and challenges in computational aerodynamics with 
emphasis on steady flow applications. Edwards and Thomas, 2 have reviewed computational methods for unsteady 
transonic flows with emphasis on applications to aeroelastic analysis and flutter prediction. Applications of CFD methods 
have become relatively commonplace in recent years, for a wide range of fluid flow problems. These applications have 
met with varying degrees of success, however, due to various deficiencies in the numerical modeling. These deficiencies 
include: an insufficient number of grid points to resolve the physics of the flow, an insufficient grid extent for external flow 
problems, too large of a time step to accurately capture the unsteadiness of the flow, and inaccurate boundary conditions, 
to name but a few. Consequently, the methods require further evaluation by performing detailed studies to assess the 
influence of each parameter in the numerical modeling on the solution. Through such studies, the accuracy and hence 
applicability of the methods may be deter min ed. 

Because of increased computer speed and memory available on current supercomputers, it is now timely to pursue 
such validation studies. Therefore, the purpose of the present paper is to report some results of such a study, performed 
using the Euler equations and the transonic small-disturbance (TSD) potential equation. Descriptions of the governing 
equations and numerical methods of solution are given first. Parameter studies for steady flow applications are presented 
next which show the effects of grid density, spatial accuracy, and far field boundary conditions. Parameter studies for 
unsteady flow applications are also presented which show the effects of reduced frequency, time step size, and the use of 
performing subiterations to minimize linearization and factorization errors in the time-accurate calculations. Comparisons 
of the TSD solutions with the Euler solutions are given to assess the applicability of the TSD potential flow method. In 
the TSD solutions, entropy and vorticity effects are accounted for to more accurately treat cases with strong shock waves. 
Finally, the results are validated by making comparisons with available experimental steady and unsteady data. 

COMPUTATIONAL PROCEDURES 

In this section, the computational procedures are described including the governing equations and the numerical method 
of solution for both the Euler equations and the transonic small- disturbance equation. 

Euler Code 

Governing Equations 

The governing equations are the time-dependent equations of ideal gas dynamics, i.e. the Euler equations, which 
express the conservation of mass, momentum, and energy for an inviscid gas. The equations are written in generalized 
coordinates and conservation form: 
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The variables £ and r; correspond to the coordinates parallel and normal to the body surface, respectively, while f 
corresponds to the spanwise direction for three dimensional flow. Q represents density, momentum, and total energy per 
unit volume. The Jacobian of the transformation, J , is defined as: 
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The equations are nondimensionalized by the freestream density Poo and soundspeed a 0 0 . 

Numerical Method 

The Euler code used in the present study is the CFL3D code developed in the Analytical Methods Branch at NASA 
Langley Research Center. 3,4 This is an implicit, finite volume, upwind-differencing code in which the spatial derivatives of 
the fluxes are split into forward and backward contributions using flux-vector splitting so that type dependent differencing 
can be utilized. The flux-vector splitting method used is that of Van Leer, 3 »” which has been recently extended for use 
on dy nami c meshes. 4 It is continuously differentiable at eigenvalue sign changes and allows shocks to be captured with 
at most two interior zones. In practice, only one zone is usually observed. 

Although flux-vector splitting is used in the current study, the code also incorporates the flux-difference approach of 
Roe. 6 The code can be used for steady and unsteady flows using either the Euler or Navier Stokes equations and can 
handle very general geometries since it allows for block structured and embedded grids. 

The starting point for the time advancement algorithm is the backward Euler time differencing scheme: 


R(Q n+l ) = o 


0) 


where R(Q n '*' 1 ) is given by: 






Q n + Jr. -Q n ~ 
w JAt w 


= 0 


(10) 


If <j> = 0, the scheme is first-order accurate in time, while if <p = 1/2, the scheme is temporally second-order accurate. 


The flux-vector splitting (FVS) method extended for dynamic meshes in Ref. 4 is used to split the fluxes into forward 
and backward contributions according to the signs of the eigenvalues of the Jacobian matrices. For example, 6^F in Eq. 
(10) at a cell centered at point * (holding the j and k indices constant) can be written: 
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The equations for F ± for a moving mesh can be found in Ref. 4. Q ± denotes state variables on cell interfaces determined 
from upwind-biased interpolations of the conserved variables. 


<?r+l/2 =«• + (JlU - + (1 + k)A+|}. (12a) 

Q .^ l /2 =<?.+! - {|[(1 - *) A + + (1 + «) A - l }. +1 ( 12 b ) 

where 

A+ = — Qi A_ — Qi — Qt-i (13) 

The parameter k € [-1,1] forms a family of difference schemes 0 : x = -1 corresponds to second-order fully-upwind 
differencing, * = 0 to Fromms scheme, and * = 1/3 to third-order upwind-biased differencing. 

When flux-limiting is desired to eliminate oscillations in shock regions, a min-mod limiter is used. Flux-limited 
interpolations are identical in form to Eqs. (12a) and (12b), except that A+ and A- are replaced with A+ and A_, 
respectively, where: 
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Equation (9) is a nonlinear equation which can be solved iteratively using a Newton linearization as 


(14) 


(15) 


^(Q e+l -Q { ) = -R(Q l ) (16) 

Here, i is a sequence of iterates so that at convergence = Q* = Q n+l . The solution of Eq. (16) at each time 
step is generally not feasible, however, since it requires the solution of a large banded matrix. The solution is therefore 
obtained with a spatially-split approximate factorization, in which the implicit operator is split into a sequence of easily 
invertible equations. The three-dimensional algorithm is implemented in three steps as: 
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For time accurate calculations, successive iterations of Eq. (17) can be used at each time step in order to minimize 
factorization and linearization errors. 


For steady state calculations, only one iteration of Eq. (17) is used at each step. In addition, the convergence rate is 
accelerated using a local time stepping procedure and the FAS multigrid scheme. 3 The implicit spatial derivatives of the 
convective and pressure terms are spatially first-order accurate, resulting in block tridiagonal inversions for each sweep. 
The algorithm is completely vectorizable since the operations in each of the spatial directions can be vectorized over the 
remaining two. 

At each time step boundary conditions are applied explicitly. In the far field, a quasi-one-dimens ional characteristic 
analysis is used to determine boundary data, while on the body, pressure and density are extrapolated from the interior 
and no flow is allowed through the surface. 
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Potential Code 
Governing Equations 

In the potential flow calculations, the flow is assumed to be governed by the general frequency, modified TSD potential 
equation which may be written in conservation law form as 


where 
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(18) 


fo = —A(j>t — B<t> x 
f\ = E<j>z + F <t> j 2 + C^y 2 
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(19) 


The coefficients A, B, and E are defined as 

A = B = 2.v4 (20) 

Several choices are available for the coefficients F, G, and H depending upon the assumptions used in deriving the TSD 
equation. The coefficients herein are defined as 

F = -±(7+l)A4 G=i(7- 3)M£ H = -h- l)Ml, (21) 

The lifting surface is modeled by imposing the following bound guy conditions: 

Flow tangency : <t> 2 ± - /j ± + ft (22a) 

Trailing wake : Ft + T z — 0 and A0* = 0 (22b) 

where A( ) represents the jump in ( ) across the wake. The flow-tangency condition is imposed along the mean plane of 
the lifting surface. In Eq. (22a) the plus and minus superscripts indicate the upper and lower surfaces of the mean plane, 
respectively. The wake is assumed to be a planar extension from the trailing edge to the downstream boundary of the 
finite-difference grid. 

Entropy effects are included in the TSD solution 7 by replacing the streamwise flux f\ (Eq. (19)) in the TSD equation 
by an alternative flux given by 
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The first term of this new flux was derived by an asymptotic expansion of the Euler equations including the effects of 
shock- generated entropy. This analysis shows that Eq. (23) is accurate to at least O(0 Z 3 ) in the expanded Euler equations. 

Vorticity effects are included in the TSD solution 7 by writing the velocity vector as the sum of potential and rotational 
components according to 


1 s 
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(25) 



In Eq. (25), the first term on the right-hand side is the gradient of a scalar potential <l> and the second term involves 
the product of the entropy s and the gradient of a Clebsch variable V. The function ^ is a measure of the stretching 
and rotating of vortex filaments associated with entropy variation. 8 For the applications of interest in the present work, 
the rotational part of the velocity vector is assumed to occur only in the region downstream of shock waves. Further 
assuming that the entropy convects with the freestream speed and that the shock curvature is negligible implies that 
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as shown in Ref. 8. These assumptions eliminate the variable $ from the formulation leaving only the entropy to be 
determined throughout the flow field. In a steady flow, entropy is constant along streamlines and changes only through 
shock waves. The entropy jump is computed along shocks using the Rankine-Hugoniot relation 
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In Eq.(28), ui is the flow speed upstream of the shock and u 9 is the shock speed. Then, for simplicity, the grid lines are 
assumed to approximate the streamlines of the flow, which is consistent with the small-disturbance approximation. The 
entropy is either convected downstream along the grid lines according to 


ds ds 
dt "** dx 


= 0 


(29) 


for unsteady applications, or is held constant along the grid lines for steady applications. 

The modified velocity vector in turn modifies the TSD equation because the streamwise disturbance speed u = 0 X is 
now given by 


U = Oz ~ 


1 3 

7(7 - l ) M£j c v 


(30) 


The new TSD equation has the same conservation law form as Eq. (18) with the new fluxes defined by simply replacing 
0 X by the modified speed given by Eq. (30). 

Numerical Method 

The potential flow code used in the present study is the CAP-TSD code developed in the Unsteady Aerodynamics 
Branch at NASA Langley Research Center. 9,10 The code uses a time-accurate approximate factorization (AF) algorithm 
similar to that described previously for the Euler equations, for the solution of the TSD equation. 11,12 The AF algorithm 
consists of a Newton linearization procedure coupled with an internal iteration technique. For unsteady flow calculations, 
the solution procedure involves two steps. First, a time linearization step is performed to determine an estimate of the 
potential field. Second, internal iterations are performed to minimize linearization and factorization errors. Specifically, 
the TSD equation (Eq. (18)) is written in general form as 


R(0 n+i ) = 0 


(31) 


where <t> n ~ rl represents the unknown potential field at time level (n-t-1). The solution to Eq. (31) is then given by the 
Newton linearization of Eq. (31) about <p l 


w ' l + ( S )../ 4=0 ,M » 

In Eq. (32), 0* is the currently available value of 0 n+1 and A 0 = 0* +1 -0*. During convergence of the iteration procedure, 
A 0 will approach zero so that the solution will be given by 0 n+1 = 0*. In general, only one or two iterations are required 
to achieve acceptable convergence. 
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The AF algorithm is formulated by first approximating the time derivative terms (fat and <p zt terms) by second-order 
accurate finite-difference formulae. The TSD equation is rewritten by substituting <t> = ^ and neglecting squares 
of derivatives of &<t> which is equivalent to applying Eq. (32) term by term. The resulting equation is then rearranged 
and the left-hand side is approximately factored into a triple product of operators yielding 




(33) 
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The equations for the spatial fluxes Fi, F 2 , F3, and the residual R are given in Ref 11. In Eq. (33) a is a relaxation 
parameter which is normally set equal to 1.0. To accelerate convergence to steady-state, the residual R may be over- 
relaxed using a > 1. Equation (33) is solved using three sweeps through the grid by sequentially applying the operators 
Lf, Lrj , and as 


£ — sweep : = -oR 

tj — sweep : Lq A<p** = A0* (35) 

£ - sweep : = A0** 

Further details of the algorithm development and solution procedure may be found in Refs. 11 and 12. 

The CAP-TSD code can be used for aerodynamic and aeroelaatic analysis of complete aircraft configurations with 
arbitrary combinations of lifting surfaces and bodies including canard, wing, tail, control surfaces, tip launchers, pylons, 
fuselage, stores, and nacelles. The code has the option of half-span modeling for symmetric cases or full-span modeling 
to allow the treatment of anti- symmetric mode shapes, fuselage yaw, or unsymmetric configurations such as an oblique 
wing or unsymmetric wing stores. 

TWO-DIMENSIONAL RESULTS 

For two-dimensional study, the NACA 0012 airfoil was selected for analysis because of the simplicity of the geometry 
and the availability of experimental unsteady data. 13 Steady-state calculations were performed for the airfoil at Moo — 0.8 
find qq = 1.25°. These conditions correspond to an AGARD test case for inviscid flow methods, 14 and have been studied 
by numerous researchers. Unsteady calculations were performed for the airfoil pitching harmonically about the quarter 
chord with an amplitude of Q\ = 2.51° and a reduced frequency of k=0.1628 (* = ) at M 00 = 0.755 and ao = 0.016°. 

These calculations are compared with the experimental unsteady data of Ref. 13. 

Steady Flow 

The effects of grid density were determined first, using the Euler code, by examining the force and moment coefficients 
as well as the pressure distribution as the grid was systematically refined. For this study, three C-type grids were 
considered. The finest grid was a 257 x 65 mesh with 257 points along the surface and wake, with 176 points on the 
airfoil surface and 65 points approximately normal to the surface. The coarsest grid was constructed by deleting every 
other mesh line from the finest grid resulting in a 129 x 33 mesh with 88 points on the airfoil. An intermediate 193 x 41 
mesh was also used which had 112 points on the airfoil. For all of the grids, the outer boundary radius was fixed at 
approximately 20 chordlengths. 

The effects of grid density on the lift, moment, and drag coefficients are shown in Fig. 1. The results are plotted as 
coefficient versus l/(total number of grid points) for three values of x, which controls the spatial accuracy of the upwind 
FVS scheme. The calculations included a point vortex representation for the airfoil to account for the induced velocities in 
the far field due to circulation (lift). As the grid density was increased the values of the lift, moment, and drag coefficients 
extrapolate to 0.3606, -0.0395, and 0.0229, respectively, for all values of x. With the third-order upwind-biased scheme 
(x s= 1/3), the lift and moment coefficients are relatively insensitive to grid density as shown in the upper part of Fig. 1. 
The solution obtained using Fromm’s scheme (x = 0) shows slightly more variation in these coefficients as the grid density 
is decreased. The second-order fully upwind scheme (x = -1) exhibits the largest change between the fine and coarse 
grid solutions. For the drag coefficient, all three schemes yield approximately the same results. With the third-order 
upwind-biased scheme, the lift and drag coefficients on the 257 x 65 mesh differ from the extrapolated values by 0.04 
percent and 1.38 percent, respectively. On the 193 x 41 mesh the values of lift and drag differ by 0.10 percent and 3.9 
percent of the extrapolated values. For both the fine and the medium grids, the moment coefficient is within 0.1 percent 
of the value extrapolated for infini te mesh density. On the 129 x 33 mesh, the lift, moment, and drag coefficients show a 
0.175 percent, 0.58 percent, and 6.8 percent difference from the extrapolated values, respectively. 

The effects of grid density on the pressure distribution are shown in the top third of Fig. 2. These results were obtained 
using the third-order scheme for the three grids previously described. The pressure distributions computed using the three 
grids are almost identical with only slight differences near the upper and lower surface shocks. 
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The effect of including the point vortex induced velocities in the far field boundary conditions on the solution is shown 
in the center of Fig. 2. The calculations were performed using the 193 x 41 grid and results were obtained both with and 
without including the point vortex. The resulting pressure distributions are very similar with small differences occurring 
only near the shocks. The effect on the force and moment coefficients, however, is more significant. For example, the lift 
coefficient decreases from the previously reported value of 0.3606, computed with the point vortex, to a value of 0.3428 
computed without it. It is noted that in Ref. 15, Euler solutions for airfoils were obtained which were virtually independent 
of the outer boundary extent when the point vortex velocities were included in the far field boundary condition. 

Finally, a comparison between Euler and TSD potential flow pressure distributions for this case is presented in the 
lower third of Fig. 2. Doth sets of results were obtained by neglecting the vortex induced velocities in the far field. The 
Euler calculation was performed using the 193 x 41 grid. The TSD calculation was performed on a Cartesian grid which 
was derived from the Euler grid by using the distribution of points along the airfoil projected to the mean plane and 
using the (approximately) vertical distribution of points at the trailing edge. As shown in Fig. 2, the two solutions are 
in very good agreement, with the main differences occuring in the regions of the shocks. The shocks are slightly sharper 
in the Euler solution, which has only one interior point in each of the shocks, in comparison with the TSD solution. 
Furthermore, this case is a very challenging one for the TSD code since it is normally regarded as being outside the range 
of applicability of TSD potential theory. The good agreement is attributable to the inclusion of entropy and vorticity 
effects as demonstrated in Ref. 7. 



2.0 r 

I 129 x 33 




Figure 1. Effects of grid density and spatial 
accuracy on the force and moment coefficients for the 
NACA 0012 airfoil at Moo = 0-9 and a 0 = 1.25° 
computed using the CFL3D code. 


Figure 2. Effects of grid density and fax field 
boundary conditions on the pressure distribution for the 
NACA 0012 airfoil at Moo = 0.8 and ao = 1.25°. 



8 


Unsteady Flow 

Results were obtained for the pitching NACA 0012 airfoil using both the Euler and potential codes. The Euler results 
were obtained using the 193 x 41 grid and the third-order upwind-biased scheme (* = 1/3). The grid for the potential 
flow calculations was derived from the Euler grid as described previously. For both Euler and TSD, the results were 
obtained using second-order time accuracy and two subiterations at each time step, unless otherwise noted. Furthermore, 
the Euler and TSD solutions did not use the point vortex far field boundary condition since the steady-state lift is nearly 
zero for this case. 

Calculated instantaneous pressure distributions at eight points in time during a cycle of motion are shown in Fig. 3 for 
comparison with each other as well as with the experimental data. In each pressure plot, the instantaneous angle of attack 
is noted. Both sets of calculations were performed using 1000 steps per cycle of motion to ensure temporal convergence. 
During the first part of the cycle, there is a shock wave on the upper surface of the airfoil and the flow over the lower 
surface is predominately subcritical. During the latter part of the cycle, the flow about the upper surface is subcritical 
and a shock forms along the lower surface. The pressure distributions indicate that the shock position oscillates over 
approximately 25 percent of the chord, and in general, that the two sets of calculated results compare well with the data 
with notable exceptions near the shocks. Furthermore, the shock from the potential solutions is generally located slightly 
downstream of the Euler shock. This result is somewhat surprising given the near perfect agreement in shock location 
for the steady case shown in Fig. 2. 





Figure 3. Comparison of instantaneous pressure distributions for the NACA 0012 
airfoil at M «, = 0.755, q 0 = 0.016°, a, = 2.51°, and k = 0.1628. 
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Effects of reduced frequency on the first harmonic components of the upper surface pressure distributions are shown in 
Fig. 4. The results are presented as real and imaginary parts of the first harmonic of the pressure coefficient, normalized 
by the amplitude of motion. Potential flow results were obtained using 1000 steps per cycle of motion for k=0.1628, 
0.3256, 0.6512, and 1.3024. Euler results for k=0.1628 were also obtained using 1000 steps per cycle while 500 steps per 
cycle were used for the other three frequencies. Comparisons between Euler and potential solutions are presented for all 
four values of k. Experimental data is available at k=0.1628 to assess the accuracy of the calculations. As shown in Fig. 
4. the two sets of calculations agree equally well across the frequency range considered and are also in good agreement at 
k=0.1628 with the experimental unsteady pressure data. 

Effects of time-step size were also investigated by obtaining results for a wide range of At using both the Euler and 
potential codes. These results are partly summarized in Fig. 5 in the form of lift and moment coefficients versus the 

instantaneous angle of attack during a cycle of motion, for the four values of k that were considered previously. Euler 
results are shown in Fig. 5(a); TSD potential results are shown in Fig. 5(b). The calculations were performed with as 
many as 1000 steps per cycle and with as few as 25. In each case, two subiterations were used at each time step. The 
Euler and TSD coefficients are generally very similar due to the similarity in the pressure distributions already shown. 
At k=0.1628, the Euler and TSD coefficients compare reasonably well with each other and with the experimental data. 
At the higher reduced frequencies, however, small differences between the two sets of results are apparennt. For example, 
at k= 1.3024 the phase lag in the Euler lift coefficient is less than that in the TSD lift coefficient. Also at k=0.6512 and 
1.3024, the TSD moment coefficients have a slightly greater magnitude in comparison with the Euler moment coefficients. 
Furthermore, the results indicate that generally much less than 1000 steps per cycle of motion were required to obtain 
temporally converged solutions depending on the reduced frequency. At the higher k-values, fewer steps per cycle of motion 
were required. At the lower frequencies, a larger number of steps were required to accurately resolve the unsteadiness 
of the flow. When too few steps per cycle were used, small differences appear in the lift and moment coefficients. The 
differences in the moment coefficients are slightly larger than the differences in the lift coefficients. These differences occur 
because of numerical oscillations which form near shock waves as shown in Fig. 6, computed using the Euler code. These 
oscillations are a further indication of an insufficient resolution of the solution in time. 




Figure 4. Effects of reduced frequency on first harmonic components of the upper 
surface pressure distributions for the NACA 0012 airfoil at Af^ = 0.755, qq = 0.016°, 
and Q\ = 2.51°. 
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Figure 5. Effects of time-step size and reduced frequency on the lift and moment 
coefficients versus instantaneous angle of attack for the NACA 0012 airfoil at 
Moo = 0.755. a 0 = 0.016°, and a x = 2.51°. 
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To examine this effect more closely, the average total variation of the pressure coefficients over a cycle of motion was 
computed using the Euler code for various numbers of steps per cycle of motion. The average total variation (TV a ) over 
a cycle was defined by 


TV a = iE„E,(|cp(« + l)-c p (i) |) (36) 

where N is the total number of steps per cycle and the summations E, and E n are over the number of grid points on 
the surface of the airfoil and over the number of time steps in the cycle, respectively. Since the total variation gives an 
indication of how much the pressure increases and decreases over the airfoil, numerical oscillations in the solution can 
be detected by a higher average total variation over a cycle of motion. Figure 7 shows the average total variation versus 
the number of steps per cycle for three values of reduced frequency. As the number of points per cycle is increased, the 
average total variation asymptotes to a particular value indicating that sufficient temporal convergence has been achieved. 
The number of steps per cycle required to adequately resolve the flow unsteadiness increases as the reduced frequency 
decreases. It is important to note however that although fewer time steps per cycle are required as the reduced frequency 
is increased, the actual time step required decreases. For example, from Fig. 7 it is seen that for a reduced frequency of 
0.1628, good results can be expected with 300 steps per cycle corresponding to a time step of approximately 0.17. For a 
reduced frequency of 1.302*4, a time step of this size yields less than 50 points in the cycle which is clearly seen in Fig. 7 
to be insufficient. This result is not surprising since the higher frequency case should require smaller time steps to resolve 
the physics of the flow. 

Finally, the effects of performing subiterations at each time step on the Euler solutions were investigated as shown in 
Fig. 8. The pressure distributions were obtained with and without subiterations using 1000 steps per cycle of motion. 
As shown in the figure, the results obtained without using subiterations are slightly different in the shock region from 
those obtained with two subiterations at each time step. If the number of steps per cycle is increased to 1700, then Fig. 9 
shows that the pressure distribution obtained without subiterations agrees closely to the previous results obtained using 
two subiterations. The influence of the subiterations is therefore to decrease the number of steps required to sufficiently 
resolve a cycle of motion. This gain, however, may be offset by the additional work required at each step. 




Figure 6. Effect of step size on the instantaneous 
pressure distribution for the XACA 0012 airfoil at 
Mto = 0.755, 0.016°, a, = 2.51°, and k =0.3256, 
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Figure 7. Average total variation versus the 
number of steps per cycle for three values of reduced 
frequency for the NACA 0012 airfoil. 


Figure 8. Effect of performing subiterations on the 
instantaneous pressure distributional «(r) - -1.25°. 



Figure 9. Instantaneous pressure distributions 
which illustrate the comparison between using a smalt 
step size with no subiterations and using a larger step 
size with subiterations. 
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THREE-DIMENSIONAL RESULTS 

For three-dimensional study, the F-5 fighter wing was selected for analysis because of the availability of both steady 
and unsteady data to assess the accuracy of the calculations. The F-5 wing has a panel aspect ratio of 1.58, a leading 
edge sweep angle of 31.9°, and a taper ratio of 0.28. The airfoil section is a modified NACA 65A004.8 airfoil which has 
a drooped nose and is symmetric aft of 40 percent chord. The calculations are compared with the experimental pressure 
data from an F-5 wing model tested by Tijdeman, et al. 16 Steady calculations were performed at Moo = 0.95 and 0.90, 
for the wing at zero degrees angle of attack; conditions which have also been studied by other researchers. Unsteady 
calculations were performed at Moo = 0.9 for the wing pitching harmonically about a line perpendicular to the root 
midchord. The amplitude of motion was selected as ai = 0.109° and the reduced frequency (based on root chord) was 
k=0.274 for comparison with the experimental pressure data of Ref. 16. It is noted that the amplitude is very small which 
may influence the accuracy of the experimental data and that the wing model had some small aeroelastic deformation 
which was not accounted for in the calculations. 


Steady Flow 

Similar to the NACA 0012 airfoil study, the effects of grid density on the steady pressure distributions of the F-5 
wing at Moo = 0-95 were determined using the Euler code. Three successively finer grids were used and the pressure 
distributions are compared with each other as well as with the experimental data. The finest grid is a 257 x 65 x 65 mesh 
which has 176 points on the wing surface in the streamwise direction, 65 points normal to the surface, and 65 points in 

the spanwise direction with 33 spanwise points on the wing. A medium grid consisted of a 193 x 33 x 41 mesh with 112 
points on the wing surface in the streamwise direction, 33 points normal to the surface, and 20 points on the wing in the 
spanwise direction. The coarsest grid, which was a 129 x 33 x 33 mesh, was constructed by deleting every other mesh line 
from the finest grid. 

Figure 10(a) show3 steady pressure distributions obtained using the 129 x 33 x 33 (coarse) and 257 x 65 x 65 (fine) 
grids at three span stations of y/s=0.18, 0.51, and 0.88, where y is the direction along the span and s is the span of the 
wing. The two sets of calculated results agree fairly well with each other except along the upper surface near the leading 
edge which indicates that the 129 x 33 x 33 grid is too coarse and that grid refinement is necessary. Solutions obtained 
using the 193 x 33 x 41 (medium) and 257 x 65 x 65 (fine) grids are compared in Fig. 10(b). At each of the three span 
stations, the calculated pressure distributions agree closely with one another at all points on the wing surface including 
the leading edge region. This indicates that both grids provide adequate grid resolution and that the results would not 
change appreciably by the addition of more grid points. Therefore, the 193 x 33 x 41 grid is adequate for obtaining 
spatially converged solutions for the F-5 wing at these conditions. Furthermore, the calculated pressures agree with the 
experimental data in regions away from the shocks. The computed shocks, however, lie aft of the experimental shock 
locations indicating that viscous effects may be important for this case. 




(a) coarse and fine mesh solutions. 



(b) medium and fine mesh solutions . 


Figure 10. Effect of grid density on the steady pressure distributions for the F-5 wing 
at Moo = 0.95 and Qo = 0.0° computed using the CFL3D code. 
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Steady pressure distributions obtained using both Euler and TSD are compared with the experimental data in Fig. 11 
for the F-5 wing at Moo = 0*90. The grid used for the Euler calculation was the 193 x 33 x 41 mesh. For the potential 
flow calculations, the grid was constructed in a similar manner as was done for the NACA 0012 airfoil. For ail three span 
stations along the wing, the agreement between the Euler and potential flow results is excellent on the lower surface of 
the wing including the leading edge suction peak. Along the upper surface of the wing, there is a mild shock wave that 
is similarly predicted by the Euler and potential codes, although the potential flow results show slightly more negative 
pressures upstream of the shock than the Euler results. The two sets of calculated results generally compare well with 
the experimental data especially along the lower surface of the wing. The data, however, does not indicate the presence 
of the mild shock wave on the upper surface at y/s=0.18 and 0.51 which is predicted by the calculations. 



Figure 11. Comparison of steady pressure distributions for the F-5 wing at 
Moq ~ 0.9 and q 0 = 0.0°. 




(a) upper surface . 



(b) lower surface. 

Figure 12. Comparison of first harmonic components of the pressure distribution for 
the F-5 wing at Moo = 0.9, a 0 = 0.0°, = 0.109°, and k = 0.274. 
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Unsteady Flow 

Results were obtained for the pitching F-5 wing using both the Euler and potential codes. The Euler calculations were 
performed using first order time accuracy, no subiterations, and slightly more than 1000 steps per cycle of motion. The 
potential flow calculations were performed using second-order time accuracy, two subiterations at each time step, and 300 
steps per cycle of motion. Figure 12 shows the real and imaginary components of the unsteady pressure distributions at 
the same three span stations as the steady results of Fig. 11. Unsteady pressures along the upper surface are shown in 
Fig. 12(a); unsteady pressures along the lower surface are shown in Fig. 12(b). On the upper surface there is a shock 
pulse in the calculated pressure distributions near 50-60 percent chord which is produced by the motion of the shock 
wave. The experimental data does not show a shock pulse in the pressures at the two inboard stations because of the 
absence of the steady shock as discussed above. On the lower surface there are positive and negative spikes in the real 
and imaginary pressure distributions, respectively, which are much more pronounced in the outboard region of the wing. 
These spikes are produced by an embedded region of supersonic flow, as indicated by the lower surface leading edge 
suction peaks shown previously in Fig. 11. In general, the two sets of calculated pressures agree well except near the 
upper surface shock pulse and in the midchord region along the lower surface. These differences may be attributed to the 
sharper shock capturing ability of the Euler code in comparison with that from the TSD code. Also, comparisons with 
the experimental data are qualitatively good for both the Euler and potential results. Further work is required to sort 
out the differences between the two sets of calculations, however, and additional correlations with experimental data are 
necessary to validate both codes. 


CONCLUDING REMARKS 

A code validation study has been presented involving the CFL3D Euler code and the CAP-TSD transonic small- 
disturbance code for both steady and unsteady applications. The purpose of the evaluation was to determine the 
accuracy and applicability of the methods by performing detailed studies to assess the influence of several parameters in 
the numerical modeling of the solution. Two-dimensional parameter studies for steady flow applications were conducted 
to determine the effects of grid density, spatial accuracy, and far field boundary conditions. The case considered was a 
NACA 0012 airfoil at a freestream Mach number of 0.8 and an angle of attack of 1.25 degrees. It was found that as the 
grid was refined, the third-order upwind biased calculation showed less sensitivity than either of the two second-order 
methods as expected. The values of the lift and drag coefficients on the finest grid (257 x 65) differ from the values 
obtained by extrapolation to infinite mesh size by 0.04 percent and 1.38 percent, respectively. The lift and drag on the 
medium mesh (193 x 41) differ by 0.10 and 3.9 percent, respectively. For both the fine and medium meshes, the moment 
coefficient is within 0.1 percent of the extrapolated value. The effect of including a point vortex representation for the 
airfoil to account for the induced velocities in the far field boundary conditions exhibited only small changes in the pressure 
distribution which when integrated over the chord, results in approximately a five percent change in the lift coefficient. 
For three-dimensional steady flow, a grid refinement study using the Euler equations was conducted on the F-5 fighter 
wing at A/oo — 0.95 and ao = 0.0° using three grids with approximately one million, twohundred fifty- thousand, and 
one-hundred forty- thousand mesh points, respectively. It was found that the finest grid and the medium grid resulted in 

essentially identical results. Comparisons between the Euler and potential flow calculations were made for both two-and 
three-dimensional cases. For the cases considered, the comparisons were very good. 

Unsteady calculations were also performed for forced sinusoidal pitching motion to determine the effects of reduced 
frequency, time step size, and subiterations performed to minimize linearization and factorization errors. It was determined 
that fewer time steps per cycle were required to resolve the flow field as the frequency of oscillation was increased. In 
addition, the influence of using subiterations at each time step is to decrease the number of steps required to sufficiently 
resolve a cycle of motion. Comparisons for two-and three-dimensional cases between Euler and potential flow calculations 
were generally good except in the immediate vicinity of the shocks where the potential flow shocks were often located 
slightly aft of the Euler shocks. This result was somewhat surprising given the excellent agreement between the Euler 
and TSD calculations for steady flow where the shocks were much stonger. 

To asses the applicability of the TSD code, comparisons with Euler calculations for both steady and unsteady flows have 
been made. The generally favorable agreement between the Euler and TSD calculations indicates that TSD calculations 
incorporating both entropy and vorticity corrections can be expected to yield results comparable to Euler calculations 
for problems similar to those considered in the current study. In addition, for TSD calculations, the grid remains 
fixed for calculating both steady and unsteady flows so that computations over complex configurations are relatively 
straightforward. Similar calculations with the Euler equations are often complicated by difficult grid generation and 
by the fact that for unsteady calculations on moving meshes, new grids must be computed at each time step. Also, 
since the solution of the three dimensional Euler equations involves five unknowns at each grid point, the computer time 
required for the Euler calculations is much higher than that required by the solution of the TSD equation. For the Euler 
code used in the present study, the computational rate for three-dimensional calculations was approximately 60 x 10~ 6 
seconds/grid point/iteration while the TSD code required only 5 x 10‘ 6 seconds/grid point/iteration. Both the Euler 
and TSD calculations were done on a CRAY-2 supercomputer at the National Aerodynamic Simulator facility located at 
NASA Ames Research Center. 
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